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Bosonic fields on rotating black hole spacetimes are subject to amplification by superradiance, 
which induces exponentially-growing instabilities (the 'black hole bomb') in two scenarios: if the 
black hole is enclosed by a mirror, or if the bosonic field has rest mass. Here we present a time-domain 
study of the scalar field on Kerr spacetime which probes ultra-long timescales up to t < 5 x W 6 M, 
to reveal the growth of the instability. We describe an highly-efficient method for evolving the field, 
based on a spectral decomposition into a coupled set of 1 + 1D equations, and an absorbing boundary 
condition inspired by the 'perfectly-matched layers' paradigm. First, we examine the mirror case to 
study how the instability timescale and mode structure depend on mirror radius. Next, we examine 
the massive-field, whose rich spectrum (revealed through Fourier analysis) generates 'beating' effects 
which disguise the instability. We show that the instability is clearly revealed by tracking the stress- 
energy of the field in the exterior spacetime. We calculate the growth rate for a range of mass 
couplings, by applying a frequency-filer to isolate individual modal contributions to the time-domain 
signal. Our results are in accord with previous frequency-domain studies which put the maximum 
growth rate at r" 1 » 1.72 x 10" 7 (GM/c 3 )" 1 for the massive scalar field on Kerr spacetime. 



I. INTRODUCTION 

An exciting era is underway in astrophysics, in which 
precision measurements of black hole masses and spins 
are becoming possible [IH1]. Naturally, there is renewed 
interest in the range of physical processes that may occur 
in the vicinity of black holes. One of the most intriguing 
processes is the so-called 'black hole bomb' effect, first 
proposed by Press and Teukolsky in 1972 [5J, in which 
the rotational energy of a Kerr black hole (BH) gen- 
erates runaway amplification in a bosonic (i.e. integer- 
spin) field. The driving mechanism is superradiance [7]: 
the stimulated emission of low-frequency radiation in co- 
rotating modes. A perturbation of frequency uj and az- 
imuthal number m, incident upon a black hole of mass M 
and angular momentum J = aM has a reflection coeffi- 
cient of greater than unity if the superradiance condition 
is met: 



UJLO < 0, 



(1) 



where uj = uj — mf2 and f2 = a/(2Mr + ) is the angu- 
lar frequency of the event horizon, with r + its radius. 
Via superradiance, a black hole may shed mass and an- 
gular momentum, but in such a way that the horizon 
area (A — 8nMr + ) actually increases. Since horizon 
area is associated with entropy [8] , superradiance may be 
thought of as a thermodynamically-favoured relaxation 
process. 

The Kerr spacetime in vacuum is (apparently) stable 
under perturbation by massless fields (see e.g. [9TU2|). 
For an instability to develop, some additional mechanism 
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is needed to confine the field in the vicinity of the black 
hole. If the BH is surrounded by a reflecting boundary 
(i.e. a 'mirror'), then superradiant flux is reflected back 
on to the hole, generating an exponentially-growing in- 
stability [S1|T3] ■ A boundary to the spacetime, as in Anti- 
de Sitter scenarios, may play a similar role |14) . A more 
'natural' astrophysical scenario arises if the bosonic field 
has a non-zero rest mass /x, in which case quasi- bound 
states [T5H2"5] may form around the BH, with Re(w) < /i. 
Bound states are localized modes which are associated 
with a minimum in the effective potential, induced by 
the gravitational attraction between the field and BH 
mass [IB]- The spectrum of bound states on the Kerr 
spacetime depends on just two dimensionless quantities: 
the BH spin a/M and the mass coupling, 



GMfx 

he 



Mil, 



(2) 



where in the following we adopt units G = c = h = 1. 

There has recently been fresh interest in a fundamen- 
tal question: is the superradiant instability in the massive 
bosonic field of any relevance to astrophysics or particle 
theory? Here, we must be mindful of a strong constraint: 
the instability only applies if M/j, < mMCl, with Mil < 
1/2, and furthermore the growth rate is exponentially 
suppressed at large m |17] , In other words, the instability 
is only physically relevant if there is a bosonic field with a 
Compton wavelength comparable to (or larger than) the 
horizon scale. Writing M/i = 7.5x 10 9 (M/M©)(/i/[leV]), 
where Mq is the solar mass, makes it clear that this con- 
dition is not satisfied by any known massive bosons in as- 
trophysical black hole scenarios. However, there are two 
more exotic scenarios which may be relevant: standard- 
model massive bosons interacting with 'primordial' black 
holes, or as-yet-undiscovered ultra-light bosons interact- 
ing with solar-mass or supermassive black holes. The lat- 
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ter scenario is motivated by compactifications of string 
theory, which suggest the existence of axions: a family 
of ultra-light bosons with masses spanning across many 
decades, possibly all the way down to the Hubble scale 
/i ~ 10 _33 eV. If extant, axions could ignite the super- 
radiant instability in stellar-mass or supermassive black 
holes, generating 'gaps' in the black hole mass-spin spec- 
trum [2HI |2"T] . It has been suggested that precision mea- 
surement of the masses and spins of supermassive black 
holes can be used to put new constraints on the existence 
of ultra-light bosons j2l)H50] . It is argued that unexpected 
gaps in the BH spectrum would be long-sought-after ob- 
servational evidence in favour of string theory (or at least, 
in favour of a beyond-Standard-Model theory requiring 
compactification of higher dimensions which leads to a 
family of ultra- light bosons). 

The majority of the literature to date has focussed on 
the instability caused by the scalar (i.e. spinless) field 
P3 [H [21 [13 IH1 EEH3S]. The massive scalar-field 
equations on Kerr spacetime are separable in the fre- 
quency domain [36] . hence a rather complete picture 
of the instability in the linear regime may be obtained 
via straightforward analy sis of ordinary differential equa- 
tions (ODEs). In Sec. II B we review the frequency- 



domain formalism and the key results in some detail. 

The aim of this work is to study the instability in 
the time domain. Time domain methods offer a gen- 
eral toolkit for studying systems of partial differential 
equations of hyperbolic character, which are well-suited 
to cases where governing equations are non-separable 
and/or non-linear. 

An interesting example of an apparently non-separable 
but linear system is the massive vector (Proca) field on 
Kerr spacetime. In the Proca-Schwarzschild case, the 
equations may be separated, by decomposing the field in 
vector spherical harmonics, leading to a single ODE de- 
scribing the odd-parity part and a pair of coupled ODEs 
for the even-parity part [371 [38]. In the Proca-Kerr case, 
impressive recent progress |291 130) has been made by 
expanding the equations in powers of a/M, applying a 
method originally developed to analyze radiation reac- 
tion in rotating neutron stars |39j . In a technical tour- 
de-force [221 [3D] , it was shown that, at first and second 
orders in a/M, the equations may once more be separated 
with vector harmonics. This has opened a window on to 
Proca-Kerr phenomenonology, which has led to a new 
constraint on the mass of the photon [29 from measure- 
ments of supermassive BHs. New time-domain methods 
[40j offer a way to move beyond the slow-rotation regime, 
and to test the applicability of the approximation. 

Non-linear effects will become important as the in- 
stability develops. An interesting question is: does the 
instability end in an explosive manner (with a 'bosen- 
ova' [2FJ), or does the system relax towards a quasi- 
stable state? A recent time-domain study [35] examined 
the evolution of an axion field with a non-linear self- 
interaction, and concluded that the system is prone to 
a sudden collapse of the axion cloud with a rapid release 



of energy. 

However, even within the linear regime, studying the 
onset of the scalar-field instability in the time domain has 
proved to be a difficult challenge, because (as anticipated 
from frequency-domain analyses [2U [33J [3T] [32]) the 
growth rate of the instability is rather slow, with an m in- 
imum e-folding time of ~ 5.8 x 10 6 M (see Sec. II B and 
Fig. [I]). To date, this timescale has proven to be beyond 
the reach of multi-dimensional simulations [351 HOI HI] , 
which (so far) have evolved the field up to t < 10 4 M or 
less. Despite this limitation, a very recent time domain 
study [3D] seems to have revealed the onset of the Proca- 
Kerr instability, confirming that it develops much more 
rapidly than the scalar-field instability. A key aim of 
this work is to construct a time-domain scheme capable 
of probing ultra-long timescales of the order of 10 6 M or 
more. We will show that ultra-long evolutions, in com- 
bination with Fourier analysis, allow us to determine the 
physically-relevant part of the bound-state spectrum to 
high accuracy. We will argue that the methods devel- 
oped here are ripe for application in the Proca-Kerr case 
(amongst others). 

Superradiant instabilities appears in a range of guises, 
for example, for BHs immersed in a magnetic field [42] . 
for higher-dimensional black holes, branes and strings 

[HGIlSim], black holes in Anti " de Sitter ( AdS ) space- 
times [HI [331 HE1 SS] and alternative spacetimes [33 and 
alternative gravitational theories [31] • The repulsive ef- 
fect of a superradiant instability on orbits around a black 
hole has also been considered [35] ■ Charged BHs gen- 
erate superradiance in charged bosonic fields which leads 
to an enhanced superradiant instability on charged rotat- 
ing spacetimes [16[l31j. On the other hand, non-rotating 
charged black holes appear to be stable under this effect 
[50] . Given the range of possible scenarios, it is possible 
that the time-domain method developed here could have 
many applications. 

In Sec. [H] we describe the basic setup and review the 
properties of bound states as revealed by frequency do- 
main analyses. In Sec. Ill we develop a time-domain 
method for evolving the massive scalar field on Kerr. We 
give a selection of numerical results in Sec. |IV| and con- 
clude with a discussion in Sec. [V] 



II. FOUNDATIONS 

In this section we recap the essential features of the 
scalar wave equation on the Kerr spacetime (II A), and 



review the properties of the bound-state spectrum in the 



frequency domain ( II B ) 



A. Spacetime, wave equation and separability 

The Kerr spacetime may be described via the Boyer- 
Lindquist coordinate system {t, r, 9, <fi}, in which the line 
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element ds 2 = g a ^dx a dx^ is 



ds 2 = - 1 - 



2Mr 



r ) 



(ft 2 



4aMrsin 2 , , , p 2 , 2 

„ dtd(j> + F —dr 2 

P 2 



-p 2 d9 2 



2Mra 2 sin 2 



A 

,2 flJ i2 



where p 2 — r 2 + a 2 cos 2 and A = r 2 — 2Mr + a 2 . The 
Klein- Gordon equation governing a scalar field $ of mass 
p is 



□$-/i 2 $ = 0, 



where 



□$ = #. a ' a = (- fl )- 1 / a ((- fl )V^$ )Q 



(3) 



(4) 



and the subscripts ';' and ',' denote covariant and partial 
differentiation, respectively. Here g = —p 4 sin 2 9 is the 
metric determinant, and g a ^ is the contravariant version 
of the metric. 

In the frequency domain, Eq. Q admits separable so- 
lutions EJ of the form 

$(t, r, 0, 0) = e- l " t e lm *i?(r)5(6»), (5) 

where the radial and angular functions satisfy 

f r + (K 2 /A A - rV) R = 0, (6) 



1 d 
sin0d0 



sm t 



d9 



A + a 2 (cu 2 - p 2 )cos 2 9 T 

sin I 



5 = 0, (7) 



with K = (r 2 + a 2 )uj — am and A = A — 2arnuj + a 2 Ld 2 . 
The solutions of the angular equation, S = Si, n (9), are 
spheroidal harmonics, and the angular eigenvalues can be 
found from a series expansion in the small-aw, ap regime 
[52] (with A — > 1(1 + 1) as auj,ap — > 0), or by solving a 
continued-fraction equation |53j . 

The radial equation can be written in standard 2nd- 
order form, by applying a simple transformation |54U55j . 
R = X/ \Jr 2 + a 2 , to obtain 



d 2 X 
dx 2 



+ U(x)X = 0, 



(8) 



with 



U(x) - (r 2 + a 2 )- 2 [K 2 - A(X + p 2 r 2 )] +p 2 -d x (3, (9) 

where f3 = rA/(r 2 + a 2 ) 2 . Here, the tortoise coordinate 
x (often denoted r*) is defined via 



dx 
dr 



a 



A 



(10) 



or, after fixing the constant of integration 
2M 



x = r + 



hi 



2M 



In 



2M 

"(11) 

In the near- horizon limit, r — > r +1 the physical solution 
satisfies the ingoing boundary condition, 



X ~ exp(— iQx), 



(12) 



where u> was defined below Eq. Q. In the far-field, the 
solution may be written as a linear superposition 



X ~ A + cxp(qx) + A exp(-qx), 



(13) 



where q = \fp? 
Re(q) > 0. 



u) 2 and we define the root so that 



B. Quasi-bound state spectrum 



The quasi-bound states are defined as the convergent 
solutions: those for which Af m = 0. The black hole is 
an open system as flux may pass into the horizon; hence 
it is natural that the quasi-bound states (like the quasi- 
normal modes [53]) have complex frequencies, 



LU = LU + iv. 



(14) 



Here, the real part ui sets the oscillation frequency of the 
mode, and the imaginary part v sets the growth (y > 0) 
or decay (y < 0) rate. The modes in the spectrum are 
labelled by integers for orbital angular momentum I, az- 
imuthal angular momentum to, overtone number n and, 
in the non-zero spin case, polarization S. If the oscillation 
frequency u) lies inside the superradiant regime (Qto < 0) 
then the bound state will grow {y > 0); if Cj = m£l then 
the mode is stationary (v = 0); otherwise the mode will 
decay (y < 0). All bound states on the non-rotating 
(Schwarzschild) spacetime are found to be decaying. We 
will confirm this picture in later sections. 

The physically-relevant bound states grow (y > 0) or 
decay (y < 0) only slowly in comparison to the Compton 
and BH light-crossing times, i.e. M|f|, \v\/fi <C 1. Foun- 
dational work by Detweiler [TH] and Zouros & Eardley 
[T7] led to simple approximations for the growth rate in 
the low- and high-Af/i regimes. These studies imply a 
maximum growth rate u max of approximately 



i((Mp) 9 (o/M - 2pr+) , Mp < 1, 
1 10 _7 e- 1 - 84M ' x , M>»1. 



(15) 



where Q is a numerical constant |68j . The maximum 
growth rate occurs in the intermediate regime Mp ~ 1/2, 
which has been explored with numerical analyses. It 
turns out that, in close analogue to the quasinormal 
modes (53j US [57], the bound states are minimal solu- 
tions of a three-term recurrence relation |23j . Finding 
minimal solutions is equivalent to finding the roots of a 
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FIG. 1: The growth rate Mv of the superradiant instability 
in the dominant mode of the scalar field (I = m = 1, n = 
0), for rapidly-rotating black holes 0.99 < a/M < 1, as a 
function of mass coupling My. The growth rates for a < M 
were found by computing the bound-state spectrum with the 
continued-fraction method of Ref. |23| . The growth rates for 
the extremal case a = M were found by direct integration of 
the radial equations. The plot shows that a maximal growth 
rate of Mv w 1.72 x 1CT 7 occurs for a w 0.997M and My w 
0.45. 



continued fraction; hence finding the spectrum numeri- 
cally is 'easy' [23]. A range of studies, by Furuhashi & 
Nambu [31], Cardoso & Yoshida [ST], Dolan [23], Rosa 
[52] and Kodama & Yoshino [25], have presented numer- 
ical results which are in accord: the maximum growth 
rate v max occurs in the fundamental bound state n = 
of the co-rotating dipole mode I = m = 1 for y near to 
|fi|. The rate v max is strongly dependent on the black 
hole spin a/M. 

Plots of the dependence of the growth rate v on a/M 
and My, were previously been presented in (e.g.) Fig. 6 in 
Ref. [33], Fig. 2 in Ref. [2T], Fig. 6 in Ref. [3T] and Fig. 14 
in Ref. [28j . Typical values of the maximum growth rate 
for a range of spin rates a < M were given in Table I 
of Ref. [33] ■ For example, the maximum rate at (e.g.) 
a = 0.99M is Mv « 1.5 x 1CT 7 , occurring at My, m 0.42. 
The extremal case a = M was considered in [32]. 

Figure [l] shows the growth rate as a function of My, 
for the dominant mode (the I = m = 1, n = bound 
state) of rapidly-rotating black holes. The plot shows 
that the maximum possible growth rate is Mf max w 
1.72 x 1CP 7 M, corresponding to an e-folding time of 
5.81 x 10 e M. Perhaps surprisingly, this maximum does 
not arise at the extremal limit (a = M), but rather at 
around a m 0.997M for Mfi 0.45. 

Under the instability, the black hole will shed mass and 
angular momentum into the dominant bound state, in the 
ratio d,J/alM = m/uj. Figure [2] illustrates the evolution 
of the black hole in the M-(J/M 2 ) plane. We make the 
simplifying assumptions that the process remains within 
the linear regime, and that the mass and angular momen- 
tum in the field are not re-absorbed by the BH as E and 



FIG. 2: Evolution of black hole parameters under superra- 
diant instability induced by a scalar field of mass fi, in the 
linear approximation. In the superradiant regime, the black 
hole loses mass M and angular momentum J into the scalar 
field (predominantly emitting into the m = Z= l,n = 
bound state). The BH parameters evolve along the coloured 
curves, from left to right. The instability ends at the marked 
points, which lie just above the (dotted) line y, = CI. The 
inset, which shows y 2 A/(2ir) as a function of J/M 2 , illus- 
trates that the horizon area A increases in this process, up 
to a maximum of A ~ iivaQ/y 2 (dotted line). [N.B. Units 
G — c = h — 1 are used]. 



J shift. The plot illustrates that, once the instability is 
in progress, the BH will migrate from left to right along 
the coloured curves, spinning-down more rapidly than it 
loses mass, until it reaches the superradiant cut-off where 
Cj = Q (~ fi). The endpoint of the instability in the lin- 
ear regime (i.e. the modes with v = 0) was the focus of 
a recent work on the extremal (a = M) case [58] , More 
realistic scenarios for the termination of the instability, 
taking into account non-linear effects, are the subject of 
recent studies f26 j 127] . [35] . 

It should be noted that the concordant view of the in- 
stability, as developed in [3TJ [53J [3_H E2] , is challenged by 
two outlier studies. First, a study of the extremal Kerr 
BH |59j predicted a maximum growth rate some four or- 
ders of magnitude faster than the consensus. However, 
it was shown in [32J that the anomalous growth rate was 
likely due to a flaw in the matching method, because the 
functional forms used in the matching procedure were not 
valid across the full range of My (the functional forms ex- 
hibited spurious poles which led to incorrect estimates). 
Second, a study of the massive scalar field in the time- 
domain in 2005 [41] reported evidence for a fast growth 
rate, Mv ~ 2 x 10" 5 (see Fig. 6 in Ref. gT]). It was 
recently argued in Ref. 40J that this finding was in fact 
due to a misinterpretation of the 'beating' effect, caused 
by interference between bound state modes of similar fre- 
quencies, which leads to a slow modulation of the enve- 
lope of the field. 
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Mirror states 



The relevant states in the case of a black hole sur- 



rounded by a mirror are defined by Eq. (12 1 and 

X(r m ) = 



(16) 



where r m is the mirror radius. The mirrored case was 
considered in detail in Ref . [T3] . 

Once more, the states are labelled by angular momen- 
tum I, m and overtone n. But, in contrast to the massive- 
field case, the number of overtones in the spectrum de- 
pends on the mirror radius r m . The number of overtones 
Ti max increases with ro, as more overtones 'fit into' the 
cavity between horizon and mirror. 

The modes are approximately evenly-spaced in fre- 
quency, with a separation that decreases as the mirror 
radius increases. For moderate-to-large r m , a simple ap- 
proximation for the oscillation frequency Cj n of the nth 
overtone is given by [T3] . 



3l+i/2,n _ (n + 1 + Z/2) 7r 



(17) 



where ji+i/2,n is the (n + l)th root of J/+i/2(0- Modes 
within the superradiant regime, 6j n < mf2, will grow with 
time, v > 0. An estimate of the growth rate is given 
in Eq. (35)-(36) of Ref. [13] . The growth rate of the 
mirrored system can be much larger than the massive- 



field case without mirror, as shown in Fig. 10 This makes 
the mirrored case rather easier to observe in time-domain 
simulations. 



III. TIME DOMAIN METHODS 

Below we develop methods for solving the scalar-field 
wave equation in the time domain. 



A. 2+1D equation 

We begin by separating out the azimuthal dependence 
of the field, 



(18) 



after introducing an alternative azimuthal coordinate tp 
via 



a 



(19) 



The new coordinate is introduced to handle the problem 
that the Boyer-Lindquist coordinate <f> is undefined as the 
horizon is approached (equivalently, to remove a term in 
the potential, a 2 m 2 , which does not vanish as r — > r + ). 
This ansatz leads to a 2+1D equation, 

{~Y?d tt + (r 2 + a 2 ) 2 d xx - AiamMrdt 

+ [2iam(r 2 + a 2 ) - 2a 2 A/r] d x -V\v = 0, (20) 



where £ 2 = [(r 2 + a 2 ) 2 - a 2 A] + a 2 A cos 2 6, and d x , d x 
denote 1st and 2nd partial derivatives, and 



V 



A 



dee 
2iam 



cot ( 



2M 
r 



1 



a 

~M~r 



+ (r 2 + a 2 cos 2 Q)ii 2 



(21) 



Eq. (20) can be solved with standard finite-difference 
methods, as in Refs. [U] ES]- Indeed, in Sec. IIIF we 



make use of a well-tested 2+ ID code to check our imple- 
mentation. 



B. Coupled 1+1D equations 



Let us now make the further decomposition, 



*= E 4>i(r)Y jm (0), 

j = \m\ 



(22) 



where Yj m (6) denotes a standard spherical harmonic 
without the azimuthal part. Now we substitute (22 1 into 



Eq. ( 20 ) and, by projecting onto the basis of harmonics, 
we recover a set of coupled 1+1D equations. The cos 2 9 



terms in Eq. ( 20 ) lead to some coupling between I modes; 



however, the coupling turns out to be rather straightfor- 
ward. We make use of the following result [ST], 



(lm\ cos 2 6\jm) 



3 jl ' 3 V 21 + 
(j, 2,0,011,0) 



2j + l 



{j,2,m,0\l,m) x 



(23) 



where 



(lm\X(0)\jm) = 27r y Yim(9)X(8)Yjm(8)d(coB0), (24) 

The numbers (j,i,m,0\l,n) are Clebsch-Gordan coeffi- 
cients. The coupling coefficient cjf is zero unless j = 1 — 2, 
I or / + 2. Hence, the even-/ and odd-/ sectors are com- 
pletely decoupled (as expected from symmetry consider- 
ations), and, within each sector, the coupling between 
/-modes of the same parity is nearest-neighbour only. 

Let us now write the equations in first-order-in-time 
form by introducing 7T; = dtipi, so that 



E 2 0) +a 2 AcI? 



tti + a 2 A (c^ +2 7r /+2 + c^_ 2 7r/_ 2 ) = 
[r 2 + a 2 ) 2 i)'{ + (2iam(r 2 + a 2 ) - 2a 2 A/r) ?/>/' 

-UamMrm -Vi- a 2 /i 2 A (c/ i ; +2 V'i+2 + 01,1-2^1-2) 

(25) 



where S 2 ^ = (r 2 + a 2 ) 2 — a 2 A and 



Vi 
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1(1 + 1) 



2AI 
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Jfr 



2iam 



p 2 {r 2 +a 2 ^)} i/n. 



(26) 
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C. Stress-energy 

The field equation ^ is obtained from an action prin- 
ciple with a Lagrangian C 



C = -\g af} d (a <S>*d p) <S> - i^ 2 $*$. 
The associated stress-energy tensor T a a is 



(27) 



(28) 



where Xr a Yp<i = ~ (X a Yp + Y a Xp). The stress-energy 
is symmetric in its indices, T a p = Tp a , and locally con- 
served, T a ^.p = 0. By contracting the stress-energy with 
a Killing vector of the spacetime we may form a conserved 
vector field. There are two linearly-independent Killing 
vectors on Kerr spacetime, and we may form vectors as- 
sociated with energy and azimuthal angular momentum 
in this way. By constructing a four-volume and apply- 
ing Gauss' theorem (as in e.g. Sec. IIB of Ref. [23]), we 
obtain two conservation laws in the following form, 



dE 
lit 



= -Ft 



dJ z 
dt 



= -T, 



(J)- 



(29) 



The 'energy' E and 'azimuthal angular momentum' J z 
are found from integrals over a timelike hypersurface. If 
we choose an 'exterior' region x > x + of a constant-i 
3-surface, we may write 



E = £dx 1 J z = 



J z dx. 



(30) 



Here £ and J z are the energy and angular momentum 
densities, given by 



(31) 



where 



2(r 2 + a 2 )£ t = A + 1) + ^(r 2 + a 2 ^)} |*,| J 
+A a 2 M 2 cr i+2 Re($r$ i+2 ) 



2 

(0) 



Aa 2 cI7) 



2Aa 2 c^ +2 Re(<i>*<i> ;+2 ) 

+(r 2 + a 2 ) 2 |$5| 2 + 2am(r 2 + a 2 )Im , (32) 

(33) 



and 



-(r 2 + a z )Ji = (S 2 0) 



=1?) mlm ( 



Acn +2 Im $i +2 <i>i + 



1+2 



-2Marm 2 \<f>,\ 2 . 



(34) 



where = r^j, 
iam<&i/(r 2 + a 2 ). 



Bt$ h *I = x $i, and $' = $' 



The fluxes are found from taking surface integrals over 
the spacelike surface at x = x+, leading to 



F [E) = (r 2 +a 2 )^Re($rl>; 



(J) 



-m(r + a") 



2 ^E 

i 



Im 



(35) 



(36) 



Let us now consider a wave of frequency to which 
is 'ingoing' in the near-horizon region, i.e. ~ 
Ai exp(— iu(t + a;)). Such a wave generates the fluxes 
Fie) ~ 2Mr+ |^4; | and Jvjj ~ (m/u))F(E)- These 
fluxes become negative in the superradiant regime < 
lu < mf2. Hence, though the wave is moving inwards, en- 
ergy and angular momentum are passing outwards into 
the exterior region. 

The choice of inner boundary x+ is somewhat arbitrary 
(typically, we take x + — — 100M), and necessitated by 
the use of a finite grid. Nevertheless, the idea of 'energy 
in the exterior spacetime' is a useful physical concept; 
we will see that, by tracking the energy and AM in the 
exterior, we may track the development of the instability. 
We note that shifting the boundary x+ merely has the 
effect of shifting the time-origin in the plots of E and 
J z , providing that x+ remains in the near-horizon region 
(x + /M«0). 



The conservation laws ( 29 1 also provide a useful code 



test, because the following quantities (representing total 
energy and angular momentum) should remain constant 
in our evolution scheme, up to numerical error: 



a 



(E) 



E- 



Jz + [ F(j)dt. (37) 



D. Finite-difference method 

We evolve the equations on a uniform grid with grid 
spacings Ax = M/n and Ai ~ Ax with n = 4 to 16. We 
use the Method of Lines, with a 4th-order Runge-Kutta 
time step, and fourth-order finite differencing on spatial 
slices, i.e. 

V'(i) ~ [-V>(i+2) + 8^(i+i) - 8^(j_i) + V'(i-2)]/(12Aa:), 
i>'(i) ~ [—0(i+2) + 16^(i+i) ~ 30^(j) 

+16^ (i -i) -^ 2) }/{12Ax 2 ) (38) 

where tj)u\ denotes ipj( x = x + + iAx) (and here i is an 
integer) . The equation set ( 25 ) resembles 



Ax = b 



(39) 



where x T = [ir k , n k+2 , ir k+4 , . . .] (where k = \m\ or 
\m\ + 1, depending on parity) and A is a real, sym- 
metric, diagonally-dominant, tridiagonal matrix. At each 
grid point, we solve (39) to find x. Fortunately, as A is 



tridiagonal, the solution can be found in a numerically- 
efficient way [6"2"| . 
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E. Perfectly matched layers 

To achieve computational efficiency, it is necessary to 
truncate the computational domain. As all physical solu- 
tions are 'ingoing' in the near-horizon region, one possi- 
bility is to impose the ingoing condition (9 t — d x ) ipi — 
on finite-difference molecules at the left boundary of 
the grid. However, this approach - imposing a condi- 
tion at an interface - runs the risk of generating high- 
frequency artifacts, which can excite numerical instabil- 
ities, and which may be mistaken for the superradiant 
instabilities themselves. In this work, we take an al- 
ternative approach: we truncate the computational do- 
main by introducing an artificial absorbing region, which 
smoothly attentuates the solution without generating re- 
flections. This is the concept of the 'perfectly matched 
layer' (PML) [53]. 

PMLs are widely used in computational electrodynam- 
ics [64] , but the underlying idea can be applied to a range 
of wave equations. For example, PMLs have previously 
been used in simulations of 'analogue' black holes in Bose- 
Einstein condensates [65} [66] . 

We may illustrate the PML approach with the simple 



example of the ID wave equation, u 



This may be 



written in fully first-order form as u — v' , v — u' . In 
the frequency domain, it has wave solutions exp(— iu)(t± 
x)). To create a PML, we may analytically-continue the 
wave solutions onto a complex- x contour and then make 
a coordinate transformation to express a; as a function of 
a real coordinate. These steps are equivalent to making 
the replacements 



d_ 
dx 



1 



d 



1 + i"f(x)/uj dx 1 



(40) 



in the frequency-domain equations. In the time-domain 
equations this leads to 



ii = v' + j(x)u, 
v = u + j(x)v, 



the solutions of which are given by 



exp ^— iu>(t ± x) ± J r y(x)da 



(41) 
(42) 



(43) 



It follows that a right-going wavepacket moving into a 
PML from the left (i.e. with y < x) will be attentuated; 
and also that a left-going wave moving into a PML from 
the right (i.e. with y > x) will also be attentuated. Note 
that here 7, the amplitude of the attentuation, is a func- 
tion of x that may be varied smoothly, over the scale 
of a few wavelengths, without generating reflections (at 
the continuum level). Thus we may 'turn on' the PML 
in a smooth manner, which greatly reduces the risk of 
numerically-generated high-frequency reflections. 

How may this idea applied to the wave equation on 
Kerr? First, we note that in the near- horizon region 
x/M <C 0, the 1+1D equations decouple, reducing to 



This equation may be turned into a simple ID wave 
equation, u — u" — by making the replacement ipi = 
exp(— ifl(t + x))u. The PML modification can then be 
applied, at the expense of introducing an auxiliary vari- 
able to play the role of v. After transforming back to the 
original variables, we are led to the set 

1p = 7T, 

7T = ip" + 2iQ(i/>' - tt) - x' ~ i^X ~ 7 (tt + i^ip) , 
X = 7 ~ X + ^VO - i^x, (45) 

(where the subscript I has been dropped for clarity). Note 
that these equations only apply in the near-horizon re- 
gion where terms ~ A are small enough to be neglected; 
this is the region where we add the Perfectly Matched 
Layer. Note also that %, the auxiliary variable, does not 
propagate outside the PML region. 



F. Validation 

We have made a number of checks of the correctness 
of our implementation. First, we tested the numerical 
convergence of the finite-different code via a ratio test, 
and found strong evidence for the expected rate of con- 
vergence. Second, we compared the results of the 1+1D 
code with those of the well- validated 2+1D code, finding 
agreement up to the anticipated level of numerical er- 
ror. Third, we tested the implementation of the PML, by 
comparing two data sets arising from evolving the same 
initial data (i) on a domain in x without a PML and suf- 
ficiently large that boundaries are not encountered, and 
(ii) on a restricted domain, x > — 200M, with a PML of 
width 10M centred around x = — 100M. Figure [3j shows 
typical results; data sets (i) and (ii) are found to be in 
excellent agreement on the right-hand side of the PML, 
as expected. Fourth, we have checked that the 'quasi- 
normal ringing' has the expected frequency and decay 
time. Fifth, we have checked that stress-energy 'con- 
stants' Cieij) [Eq. (37)] are preserved in our runs, up 
to numerical error (see Fig. 11). Sixth, we have checked 
that the results presented here are robust under adjust- 
ments of 'internal' parameters (such as the width and 
amplitude of the PML; the maximum number of I modes 
used; and the grid resolution). Finally, the next sections 
show excellent agreement with frequency-domain results, 
which we take to be further evidence in favour of a valid 
implementation. 



ipi - V" + 2imQipi - Hm£h\)\ = 0. 



(44) 



IV. RESULTS 

In this section we present a selection of numerical re- 
sults, for two cases in which an instability is present: the 
massless scalar field with a mirror (Sec. |IV A ), and the 
massive scalar field without mirror (Sec. IV B). In each 
case, we seek to stimulate the system with some broad- 
band initial data - typically a Gaussian in x - and study 
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FIG. 3: Validation of the Perfectly-Matched Layer method. 
Plot (a) compares snapshots taken at t — 140M, from two 
evolutions of the same initial data (a time-symmetric Gaus- 
sian of width o x = 5M centred around xo = 10M in mode 
I = m = 1) for a BH of spin a — 0.99A/, using (i) [red, 
solid] no PML, and (ii) [blue, dashed] a PML centred around 
x — — 100M of amplitude 7(1) (with 7/IO shown as dashed 
magenta line). Outside the PML region, (i) and (ii) are in 
good agreement. Inside the PML region, incident waves are 
smoothly attenuated without reflection. Plot (b) shows the 
same data on a logarithmic scale. Plot (c), showing data ex- 
tracted at x — x\ = as a function of time, illustrates that, 
to the right of the PML, data sets (i) and (ii) are virtually 
indistinguishable. 



perturbation. Figure [4] compares the radial profile of the 
response of a Schwarzschild BH (a — 0) with that of 
a rapidly-rotating Kerr BH (a = 0.99A/), over several 
decades of simulation time t/M = 10, 100, . . . 10 5 . Here 
the initial data is a time-symmetric Gaussian, 

ij l=1 =exp(-(x-x )/[2al)), ^ = = 4 (46) 

with xq = 10M and a x = 2M. At early times (e.g. t = 
10M), the response of the two BHs is very similar. At 
intermediate times, t = 100A/-1000M, some part of the 
initial perturbationhas been absorbed through the hori- 
zon, and some remains in the exterior; the responses 
of Schwarzschild and Kerr are now quite different, but 
the field amplitudes remain comparable. At late times, 
t = 10 4 M and t = 10 5 Af we see that the scalar field 
on the Schwarzschild BH has essentially disappeared 
through the horizon. On the other hand, the scalar field 
on Kerr has experienced significant amplification (note 
vertical scale). 

These conclusions are supported by Fig. [5] which shows 
the radial profile of the energy density 8 [Eq. (32)] com- 
puted from the field shown in Fig. |4| The late time plots 
suggest that there is more than one growing mode in the 
Kerr case, and that the mode with a peak closest to the 
mirror surface is emerging to dominate the signal. 

Let us now consider the field measured at a fixed point 
x = x\ as a function of time. Figures [6] and [7] (log scale) 
show typical time series data for Schwarzschild and Kerr 
BHs. In the Schwarzchild case, the field decays exponen- 
tially, as the mirror reflects flux back onto the absorbing 
horizon. In general, the decay rate increases as the mir- 
ror is moved closer to the BH. In the Kerr case, the field 
is long-lived. Figure [6] shows that there is typically a 
'beating' effect in the signal, generated by interference 
between modes of different frequencies. This leads to a 
modulation of the signal which can disguise the under- 
lying growth and even lead to incorrect conclusions (as 
argued in Ref. [40J. On the other hand, Fig. [7| shows 
that the total energy in the exterior E (Sec. Ill C ) gives 
a clear view of the underlying growth as: (i) E does not 
oscillate sinusoidally, and (ii) the variation of E depends 
only on the energy flux through the inner boundary, Ft E ^ 
[Eq. ( 35 )] , which is only marginally influenced by beating 
between modes. 



the features of the response. In particular, we attempt 
to compute the spectrum of superradiant bound states 
which dominate the response at late times. We will focus 
on the m = 1 mode, on which superradiance is expected 
to have the strongest effect. 



1. Mirror spectrum 

The spectrum of the response is revealed by taking a 
Fast Fourier Transform of the time series data for the 
field (extracted from the grid at fixed x — X\). Let us 
introduce the (complex) Fourier amplitude, defined by 



A. Scalar field with mirror 

Let us begin by considering a BH surrounded by a 
mirror at r m = 20M, subject to some smooth initial 



N-l 



fl(vj) = V^( i fe) ex P(-^ i i fe ), 



(47) 



fe=0 



where — fcAi, ujj — jAu, and N is the number of data 
points, At = t m!iX /(N — 1) and Aw = 27r/i max , and the 
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FIG. Ax Time evolution of field on black hole background with 'mirror' at r m = 20M. These plots show snapshots of ipi 
[Eq. ( 22 1] for the dipole mode I = m= 1 of the scalar field, as a function of tortoise coordinate x/M, at t = (top left), 10M, 
100A/, 10 3 M, 10 4 M and 10 5 Af (bottom right), for the Schwarzschild (a = 0, solid line) and fast-rotating Kerr (a = 0.99M, 
broken lines) BHs. Here the initial data is a time-symmetric Gaussian ( |46[ ) of width a x = 2M centred about xq = 10M. By 
late times, the field has been absorbed by the Schwarzschild BH, and amplified by the Kerr BH. 



power spectrum given by 



P l {w) = \f l {u)Y 



(48) 



Note that the frequency resolution Aw is inversely pro- 
portional to the simulation time i max ; hence longer runs 
yield a more refined view of the spectrum. 

Figure [8] shows typical power spectra, for mirror radii 
r m = 15M and r m — 25M. The sharp peaks indicate the 
presence of modes in the spectrum. The plot makes it 
clear that the modes within the superradiant regime < 
Moj < mil have the largest amplitude (highest peaks), 
as expected. 

The spectrum has all the features anticipated in the 
frequency-domain study of the mirrored case, Ref. |13j . 
The peaks are approximately evenly-spaced in frequency, 
with a separation that decreases as the mirror radius in- 
creases, as predicted by Eq. (17). 



The width of a peak in the power spectrum provides 
information on the decay/growth rate of the correspond- 
ing mode. Consider a monochromatic mode with time 



dependence cxp (— i(uit - 
power spectrum is 



iv)t). Its contribution to the 



P(w) sa 1/ [{uj-Lu f + , 



(49) 



Hence, fitting a quadratic to l/P{oS) in the vicinity of a 
minimum w should give an estimate of the decay/growth 
rate whose sign may be found by examining /(w)). De- 
spite its promise, we found that this method was not suffi- 
ciently accurate to numerically estimate the rates of long- 
lived modes, and so we pursue an alternative method 
below. 



2. Growth rates from frequency-filter 

To study the individual modes in the time-series data, 
we devised a simple 'frequency-filter' method. First, we 
identified the modes in the power-spectrum (e.g. Fig. [8]) 
using a peak-finding algorithm, and used interpolation to 
obtain a better estimate for the real frequencies u> n . Next, 
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FIG. 5: Time evolution of energy density £ on black hole spacetime with 'mirror' a t r m = 20M. These plots show snapshots 
of the energy density obtained from the field shown in Fig. [4] (see caption) via Eq. ( 32 1 . The last plot (note the scale) shows 
that the energy of the field on the rotating BH (broken line) has increased by a factor of ~ 100 over a timescale of ~ 10 5 M via 
the superradiant mechanism; whereas the field on the non-rotating spacetime (a = 0) has drained away through the horizon. 



for each mode, we (i) multiplied the Fourier transform by 
a narrow filter function centred on the peak, i.e. f(u>) —> 
f(u))F(u) — Li n ) with 



F(x) = cxp(— x 4 /w 4 ) 



(50) 



and w — fcAw with n ~ 5 (typically) and Aw defined 
above Eq. ( |48[ ) ; (ii) reconstructed the modal contribu- 
tion by taking the inverse Fourier transform of the filtered 
data; (iii) computed the square magnitude of the recon- 
structed field as a function of time; and (iv) computed 
the growth index using linear regression on the logarithm 
of the field magnitude. 

These steps are illustrated in Fig. [9j which shows the 
square magnitude of the reconstructed field modes for 
the first few overtones. It shows that the reconstructed 
field modes at early and late times depends somewhat on 
the shape of the filter; but that the field around t max /2 
is a faithful reconstruction of the modal contribution. 

With the frequency-filter method, we may obtain accu- 
rate estimates for the frequency and growth/decay rate 
of the subdominant modes that are 'hidden' in the time- 



domain data. Figure [TU] shows estimates of the growth 
rate of the instability, as a function of mirror radius, 
for the first few overtones. These estimates are com- 
pared with the growth rates obtained in Ref. [131 via a 
frequency-domain analysis. The agreement is found to 
be excellent. It is gratifying that many overtones can be 
extracted from the time domain data, even in cases when 
the late-time signal is dominated by the fastest-growing 
mode. 



B. Massive scalar field 

In this section, we consider the case of a massive scalar 
field in vacuum. Non-zero mass leads to a spectrum of 
quasi-bound states with Re(w) < fx, as summarized in 



sec, una 



Figure 11 shows the evolution of the energy and angu- 
lar momentum of the massive scalar field in the exterior 
spacetime, for a typical simulation of the m = 1 mode of a 
scalar field with M/x = 0.42 and Gaussian initial data, on 
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t/M 



FIG. 6: Time series data for the evolution of the field on a BH 
with mirror. The plot shows the (real part of) the field shown 
in Fig. [4] and [5] extracted at r = WM, as a function of time 
t/M. On the non-rotating spacetime (a = 0, blue) the field 
decays, whereas on the rotating spacetime (a = 0.99M, red) 
it is long-lived. The latter case shows evidence of 'beating', 
due to interference between modes of different frequencies, 
and growth (see Fig. [7] . 
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FIG. 7: These plots show the time evolution of the field as 
in Fig. [fi] on a logarithmic scale, up to t = I0 5 Af. The up- 
per plot shows that the envelope of the Kerr (Schwarzschild) 
field grows (decays) exponentially. The lower plot shows that 
the energy in the exterior (see Sec. Ill C I also grows (decays) 
exponentially without beating effects, and with twice the ex- 
ponent, as expected. 



a Kerr BH with a = 0.99M. At 'early' times t < 10 4 M, 
which have been probed by other studies [23 SOI EI], 
we see that the energy in the exterior initially drops off, 
as non-superradiant flux is absorbed by the BH, before 
approaching an apparent 'stationary' equilibrium. How- 
ever, 'late' time evolutions t 10 4 M, reveal that the 
solution is not truly stationary, as the slow exponential 
growth in the bound states begins to dominate. 

To check the validity of our results at late times, we 
have computed the energy and AM 'constants', defined 
in Eq. (37), from the numerical data. Figure 
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c) shows 

that there is a slow 'drift' in the numerical values of 
the constants, but that the drift is strongly resolution- 
dependent, as expected. We also see that the growth 
rate is resolution-dependent, but that it converges to a 
physical value as n is increased. 



1. Bound state spectrum 

By extracting the field at a fixed point x = x%, and 
taking the Fourier transform, we may obtain the power 
spectrum Pi(ui) of the massive field (see Sec. IV A 1 ). Fig- 



ure [12] shows a typical spectrum in the case a = 0.99A/, 
M/j, = 0.42. The fundamental mode and the first four 
overtones, n = . . . 4, are apparent as sharp peaks in the 
spectrum. With ultra-long timescalc simulations, we gain 
sufficient frequency resolution to determine the bound- 
state frequencies with high accuracy. Figure [12] shows 
a comparison between the frequencies found by Fourier 
analysis of the time domain signal, and the frequencies 
found from frequency-domain analyses (e.g. Sec. II B and 
Ref. [13]). 

By applying the frequency-filter method described in 
Sec. |IV A~T] we may estimate the growth rate of the super- 
radiant bound states. A key advantage of the method is 
that it allows us to separately analyse the many overtones 
that are in the spectrum, even though the time-domain 
signal may be dominated by a single mode. Further- 
more, we find that we may extract good estimates for 
the growth rate of lowest overtones n = 0, 1 from runs of 
i max = 10 5 M which are relatively 'short' in comparison 
to the e-folding time (Mu > 5.8 x 10 6 M) of the insta- 
bility. Figure [13] shows some typical results, comparing 
the growth rate extracted from time-domain simulations, 
with the much more accurate values found by frequency- 
domain analysis |23j . The agreement is good, reinforcing 
once more the concordant view of bound states described 
in Sec. iHBl 



V. DISCUSSION AND CONCLUSION 

In this paper we have described a 'coupled 1+1D' time- 
domain scheme for evolving the scalar field on Kerr space- 
time. The scheme enabled us to study the evolution of 
the scalar field, within the linear regime, over ultra-long 
timescales t ~ 10 6 M, which represents a factor of 100 
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FIG. 8: Power spectrum of the scalar field on BH with mirror. The plot shows the power spectrum (48 1 as a function of frequency 
Mui, computed via a Fourier transform of the time-domain data (Fig.[7|, for two mirror radii, r m = 15M and r m = 25M. Modes 
inside the superr adia nt regime, < u < mfl, are subject to exponential growth. The modes are approximately evenly-spaced 
in frequency [Eq. [l7] . The plot illustrates that the effect of increasing the mirror radius r m is to bring more overtones within 
the superradiant window. 




domain analyses (see Fig. 10 and Fig. 13). Thus, our 
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FIG. 9: The frequency-filter method. Each solid line in the 
plot corresponds to the (normalised) square amplitude of a 
particular mode extracted from the time-domain data. Each 
mode is isolated by applying a narrow frequency filter to the 
Fourier transform of the time domain signal, and then apply- 
ing the inverse transformation (see text). The dotted lines 
represent the best fits. The plot shows the overtones n — 
to 3 for r m = 32M and a = 0.99M, for which case n = 3 is 
the most rapidly-growing. 



to 1000 improvement over previous time-domain studies 
[23 |40l |41] . The efficiency of the scheme allowed us to 
track the massive-field superradiant instability for a sig- 
nificant fraction of an e-folding time. We have shown 
that accurate values for the bound state frequencies and 
growth rates may be extracted directly from time-domain 
data via Fourier analysis. We found that these values 
were in excellent agreement with the results of frequency- 



technique offers a promising way forward for tackling sce- 
narios where a frequency-domain analysis, reliant on the 
separability of the underlying equations, is not feasible. 
The method would seem to be ripe for application to 
the Proca-Kerr system [29, 30, 37], recently examined by 
Witek et al. [3D] in the time domain for the first time. 

The efficiency of the method is primarily due to the di- 
mensional reduction procedure, i.e., working in a 1+1 D 
rather than 3+1D domain. In other words, the key step 
in our method is the decomposition of the field into a su- 
perposition of harmonics. Even though this does not lead 
on to a full decoupling, it turns out that the coupling be- 
tween I modes is very simple. In the scalar-field case, the 
coupling is only between the nearest-neighbour modes of 
the same parity. In the Teukolsky equation (s ^ 0), we 
anticipate the coupling would also be simple, but instead 
between nearest- and next-to-nearest neighbours [61] . In 
the higher-spin massive-field case (e.g. Proca) the cou- 
pling may be substantially more complicated, although 
the slow-rotation expansion of 29 , 30 suggests it may be 
tractable. 

Our 1+1D scheme is similar in some regards to that 
developed in Ref. [1)7], to look at the late-time tails of 
massless fields on Kerr spacetime. In Ref. [57], the coef- 
ficient of dtt^ , which depends on cos 2 9, appears in the 
denominator on the right-hand side of the equation. The 
denominator is converted to an infinite series in the nu- 
merator, which leads to coupling between all I modes 
after harmonic decomposition (see Sec. 2.2 in 67J). For- 
tunately the couplings die away rapidly with increasing 
I and the method is useful in practice. An advantage 
of the approach of Ref. |67j is that a matrix-inversion 
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FIG. 10: The growth rate of the scalar field for rotating BH with mirror. The plot shows the growth rate of the first few 
exponentially-growing modes (oc exp(i^)), as a function of mirror radius r m . The points show the gro wth rates calculated 
from time-domain data, using runs up to t = 10 M, by applying the frequency-filter method (Sec. IV A 2 1. The lines show the 
growth rates found in the frequency-domain analysis in Ref. |13j (with thanks to the authors). 



step (Sec. HID) is not required at each time step. This 
may offer an alternative route to tackling the Proca-Kerr 
case, where coupling terms may be substantially more 
complicated. 

Our method also makes use of a technique from 
computational electrodynamics: the Perfectly Matched 
Layer [531 IM] ■ We have shown that the PML offers a 
numerically-robust way of implementing the absorbing 
boundary condition at the event horizon. It seems that 
this technique may be readily incorporated into a range 
of linear wave equations, and yet to date it not been 
used widely in black hole perturbation theory. The PML 
also offers distinct advantages in multi-dimensional sce- 
narios, where stable boundary conditions are more diffi- 
cult to implement. It would be interesting to compare 
the efficacy of the PML with alternative methods, such 
as horizon-penetrating coordinate systems with excision. 

We have run our time-domain code for a range of mass 
couplings, 0.1 < 1, and BH spin rates < a/M < 0.99, 
for low azimuthal numbers m. In every case, we have 
found that the growth of the field (if present) is due 
to growth in a superposition of bound states which lie 
within the superradiant regime. Furthermore, these 
bound states are found to have the expected frequen- 
cies and growth rates from the 'concordant' frequency- 



domain analyses (see Sec. II B). We have found no evi 



dence in favour of the more rapid growth rates reported 
in [4l] and [59]. 

Unabated, an exponentially-growing instability will in- 
flate the magnitude of the scalar field until non-linear ef- 
fects can no longer be ignored. Non-linear effects arise 
from (at least) two sources: the self-interaction of the 
field, and the back-reaction of the field upon the space- 
time. In the context of axions, the former effect may be 
modelled by replacing /i 2 $ with /^ 2 sin$ in the scalar- 
field equation, introducing a non-linear term at leading 
order /z 2 $ 3 . A recent pioneering time-domain study |35j 



looked at the effect of non-linear coupling. It concluded 
that, whereas small self-couplings will increase the rate 
of energy extraction, large self-couplings will cause a col- 
lapse of the axion cloud. This suggests that the end- 
point of the instability is an explosive phenomenon |26| . 
As yet, the transition from small to large couplings has 
not been tracked via a time-domain simulation, due to 
strong limitations on run-time imposed by the computa- 
tional expense of 3+1D codes. It may be possible that 
a dimcnsionally-reduced code could allow longer run- 
times; however, the non-linear terms would introduce 
couplings between m and I modes, which would substan- 
tially increase the complexity of implementation. A per- 
haps more tractable problem would be to use a time- 
domain code to study the generation of gravitational 
waves (i.e. linearized metric perturbations on the Kerr 
background) induced by the stress-energy of the scalar 
field. 

A further prospect is to apply dimensional reduction 
to study the Proca field over long timescales within the 
linear regime. Here a 2 + ID code, exploiting a simple 
separation into m-modes, would seem suitable. A 2 + ID 
code would allow one to probe beyond the timescales 
achieved in Ref. [30] (£ max = 3 x 10 3 M). The resolution 
of the power spectrum scales with run-time, and so, with 
long-timescale datasets, one would be able to isolate the 
modes by applying a narrow filter to peaks in the power 
spectrum (as in Sec. IV A2) to determine the details of 



the bound states (i.e. frequencies and growth rates) to 
high precision. This is surely an interesting challenge. 

In summary, time-domain simulations offer a way to 
examine the onset, development and termination of black 
holes superradiant instabilities which may be ignited in 
nature by ultra-light bosonic fields. It is important to 
understand these instabilities fully, if we are to use new 
observations of astrophysical black hole parameters to 
constrain the boson mass spectrum [26-30, 35J. Further- 
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Time-domain evolution of massive scalar field on Kerr: 'early' phase 
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Time-domain evolution of massive scalar field on Kerr: 'late' phase 
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FIG. 11: Time-domain evolution of a massive scalar field {M\i = 0.42) on rapidly-rotating BH spacetime (a — 0.99M) with 
time-symmetric Gaussian initial data, with x — 20M and a x = 10M [see Eq. ( 46 1] . The plots show the evolution of the energy 



E and azimuthal angular momentum (AAM) J z in the exterior spacetime ( 30 1, and the 'conserved' quantities ( 37 1. In the initial 
phase, t < 10 4 M shown in plot (a), E decreases and J z increases as non-superradiant flux, particularly th e counter-rotating 
part of the field with uj < 0, is absorbed by the BH. The 'total' energy and AAM, C^e) an d C(j) [Eq. ( 37 1] , are numerically 



conserved to a good approximation. In the 'late' phase, t > xl0 4 M shown in plot (b), the exterior energy and AAM increase, 
as growth in the superradiant bound states begins to dominate. Plot (c) shows that there is some numerical violation of the 
constraints which diminishes as resolution improves, and which is consistent with the resolution-dependence of the curves at 
late times seen in plot (b). 
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FIG. 12: Power spectrum of massive scalar field (Mn = 0.42) on fast-rotating BH spacetime (a = 0.99M) in the co-rotating 
dipole I = m — 1. The spectrum is computed from time-domain response to Gaussian initial data, up to t — 6 x 10 M. The 
fundamental mode and the first four overtones are visible. The frequencies of the first bound states, n = . . .4 for I = m = 1, 
computed by frequency-domain analysis, are shown at uj/fi w 0.9733, 0.9883, 0.9936, 0.9960 and 0.9980. 
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FIG. 13: Instability growth rates of massive scalar field on rapidly-rotating Kerr spacetime (a — 0.99M). The lines show 
the growth rates of the fundamental and first-overtone bound states of the co-rotating dipole I — m — 1, established via the 
frequ ency-do main method [23]. The points show estimates of the growth rates found by applying the frequency-filter method 
(Sec. IVA2I to data extracted from time-domain runs (with t max = 10 J Af and resolution dx = M/6). The observed agreement 



implies that the bound-state spectrum can be accurately determined from 'ultra-long' time-domain datasets in cases, such as 
the Proca field on Kerr spacetime, where a lack of separability impedes a frequency-domain analysis. 



more, the study of field configurations which are long- 
lived or quasi-stable is of foundational interest, because 
such configurations would seem to challenge to the spirit 
of the 'no-hair' conjecture. We have shown here that the 
scalar-field instability has been well understood within 
the linear regime. There remains much scope for pursu- 



ing the instability into the non-linear regime, and in the 
context of massive vector fields. 
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